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Results from twin control simulations of the preindustrial C0 2 gas exchange (natural flux of C0 2 ) between 
the ocean and the atmosphere are presented here using the NASA-GISS climate model, in which the same 
atmospheric component (modelE2) is coupled to two different ocean models, the Russell ocean model 
and HYCOM. Both incarnations of the GISS climate model are also coupled to the same ocean biogeo- 
chemistry module (NOBM) which estimates prognostic distributions for biotic and abiotic fields that 
influence the air-sea flux of C0 2 . Model intercomparison is carried out at equilibrium conditions and 
model differences are contrasted with biases from present day climatologies. Although the models agree 
on the spatial patterns of the air-sea flux of C0 2 , they disagree on the strength of the North Atlantic and 
Southern Ocean sinks mainly because of kinematic (winds) and chemistry (pC0 2 ) differences rather than 
thermodynamic (SST) ones. Biology/chemistry dissimilarities in the models stem from the different 
parameterizations of advective and diffusive processes, such as overturning, mixing and horizontal tracer 
advection and to a lesser degree from parameterizations of biogeochemical processes such as gravita- 
tional settling and sinking. The global meridional overturning circulation illustrates much of the different 
behavior of the biological pump in the two models, together with differences in mixed layer depth which 
are responsible for different SST, DIC and nutrient distributions in the two models and consequently dif- 
ferent atmospheric feedbacks (in the wind, net heat and freshwater fluxes into the ocean). 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

Over the 20th century the ocean has taken up about 1.5- 
2.0PgCyr _1 , i.e. about 25-30% of all the manmade emissions of 
C0 2 from fossil-fuel burning, cement manufacturing and land use 
change (Khatiwala et al., 2009; LeQuere et al., 2009; Sabine et al., 
2004). The geographic pattern of the oceanic sink is characterized 
by uptake at higher latitudes and outgassing in the tropics (Takah- 
ashi et al., 2009; Gruber et al., 2009). The ocean-atmosphere C0 2 
exchanges, which act on timescales of hours to days, and the 
deep-ocean export of carbon, which acts on timescales of days to 
years and longer, together with land-atmosphere exchanges, set 
the atmospheric C0 2 concentration and its radiative impact on 
Earth’s climate (Friedlingstein et al., 2006). At the same time, ocean 
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circulation changes due to variations of the physical forcing of the 
ocean, such as winds and air-sea freshwater fluxes, may affect the 
geographical patterns of deep water formation leading to large 
changes in the C0 2 uptake patterns (Marinov et al., 2008b). 

Several studies have emphasized the uncertainty in the esti- 
mates of the magnitude and variability of the air-sea flux of C0 2 
in the Earth’s climate system (LeQuere et al., 2009; Gruber et al., 
2009; Takahashi et al., 2009; Friedlingstein et al., 2006). These 
uncertainties mainly arise from incomplete observational coverage 
but also from insufficient model complexity, i.e. not realistic pro- 
cess simulation or lack of components, as well as model errors. 
For example, the interannual variability and the long-term trends 
in the flux of C0 2 in the atmosphere, either the natural trends or 
the anthropogenic trends, are affected primarily by changes in 
sea surface temperature (SST), surface wind and ocean pC0 2 distri- 
bution. Ocean biogeochemistry (the biological or soft-tissue pump) 
affects dissolved inorganic carbon (DIC) which in turn determines 
pC0 2 . Doney et al. (2004) showed that model biases can be attrib- 
uted to changes in surface forcing, subgrid-scale parameterizations 




A Romanou et al./ Ocean Modelling 66 (2013) 26-44 


27 


and model complexity, all of which affect model-ocean transport 
and dynamics and eventually degrade ocean tracer and carbon cy- 
cle variable distributions (e.g. air-sea C0 2 flux, anthropogenic C0 2 
uptake and export production). Marinov et al. (2008a) illustrated 
how deep ocean ventilation determines the pathways that atmo- 
spheric C0 2 follows after uptake and how these pathways are af- 
fected by the biological pump. 

Uncertainties become even more apparent in a coupled atmo- 
sphere-ocean-carbon model setting in which errors in the ocean/ 
atmospheric background models can be amplified or diminished 
by the coupling and hence mask the real climate/carbon-cycle 
feedbacks in such models (Cox et al., 2000). 

Assessing the model carbon cycle and the uncertainties associ- 
ated with the representation of physical and biogeochemical pro- 
cesses is crucial in improving the models and subsequently the 
estimates of the change in air-sea carbon fluxes. The physical 
model describes the large-scale vertical exchanges (ventilation, 
convection and water mass formation, subduction) and the lateral 
transports (both mean and eddy induced) that affect the biogeo- 
chemistry response and ultimately the gas exchange flux to/from 
the atmosphere (Gnanadesikan et al., 2002; Schneider et al., 
2006; Ito et al., 2010). These processes affect the distributions of 
preformed DIC in the ocean (Marinov et al., 2008a) as well as 
the upward transport of nutrients into the euphotic zone needed 
to sustain the marine biological production. Downward transport 
(export) of particulate organic matter from the upper water col- 
umn to depth also changes the vertical gradient of DIC in the 
water column (Volk and Hoffert, 1985). Deep carbon export flux 
determines the efficiency with which carbon is transported to 
depth away from further contact with the atmosphere on long 
time scales (Henson et al., 2012). Variations in the export lead to 
changes in the air-sea C0 2 exchanges and the partitioning of car- 
bon between the oceanic and atmospheric reservoirs (Schneider 
et al., 2006). 

The purpose of this paper is to examine the sensitivity of the 
carbon uptake by the ocean to parameterizations that relate to dif- 
ferent vertical coordinate systems, thereby assessing the uncer- 
tainty in the carbon sources/sinks distribution based on physical 
model differences. This is particularly interesting, albeit complex, 
in the context of coupled atmosphere-ocean-ice models where 
feedbacks may lead to error compensation or amplification in the 
determination of the C0 2 fluxes. Specifically, we examine the 
divergence in model solutions during spin up and at equilibrium 
of two carbon cycle climate models that share the same atmo- 
spheric, ice and ocean biogeochemistry components while differ 
in the ocean model component. This assessment will set the stage 
for a better determination of the biases in the transient, present 
day climate, and future climate simulations. 

The paper is arranged as follows: Section 2 provides a descrip- 
tion of the two GISS climate models with emphasis placed on the 
oceanic components and parameterizations therein, the biogeo- 
chemical model, the gas exchange and the radiative 3-way cou- 
pling between the ocean, the atmosphere and the 
biogeochemistry. Section 3 provides an overview of the external 
datasets used in this study. These are either observational or syn- 
thetic datasets used for forcing, initialization and model validation. 
In the fourth section we present the model results; firstly focusing 
on the physical model differences with present-day-observed cli- 
matologies in order to assess the realism of the simulations, the ex- 
tent of the inter-model differences, as well as elucidate the reasons 
that such differences exist. The second part of the Results section 
describes and evaluates the carbon cycle differences between the 
models. This part is organized by the three major pathways carbon 
is distributed in the oceanic system, the solubility, biology and car- 
bonate pumps. Finally, Section 5 summarizes the conclusions of 
this study. 


2. Model description 

Two incarnations of the GISS climate model exist coupled to 
prognostic ocean carbon cycle, both using the same atmospheric, 
land and ice components but different ocean models. Results from 
these runs have been submitted to the Coupled Model Intercom- 
parison Project Phase 5 (CMIP5; http://cmip-pcmdi.llnl.gov/ 
cmip5/) as GISS-E2-R-CC and GISS-E2-H-CC which are the carbon 
cycle (CC) versions of modelE2 coupled to the Russell (R) or Hycom 
(H) oceans. In this paper, we will be using the abbreviations GISSER 
and GISSEH to refer to GISS-E2-R-CC and GISS-E2-H-CC respec- 
tively. A description of each component is provided below, with 
the emphasis placed mainly on the carbon model and the differ- 
ences between the two ocean models. 

2.2. The atmospheric model, modelE2 

ModelE2 is an updated version of the GISS atmospheric model 
(modelE) that was included in the Climate Model Intercomparison 
Project (CMIP3) (and IPCC 4th Assessment Report) simulations 
(Schmidt et al., 2006; Hansen et al., 2007) and which is currently 
used in the CMIP5 experiments (Schmidt et al., 2012). 

Improvements include (a) a finer horizontal resolution, at 2 by 
2.5° latitude/longitude and 40 levels in the vertical with the top 
of the atmosphere at 0.1 mb, (b) refinements in the cloud micro- 
physics, (b) a new turbulence scheme with no dry convection, (c) 
a new gravity wave drag parameterization, (d) new glacial-melt 
and melt-pond fraction and (e) updated radiative transfer, aerosol 
and ozone input datasets. Model time steps are about 4 min for the 
core dynamics, 30 min for the model physics and 2.5 h for radiative 
calculations in order to reduce computational costs. The model em- 
ploys a quadratic upstream differencing scheme for the advection 
of tracers such as C0 2 that is very accurate and non-diffusive 
(Prather, 1986). 

2.2. The ocean models 

ModelE2 is coupled to two different coarse-mesh ocean models, 
either to the Russell Ocean model (GISSER) or to HYCOM (GISSEH). 
The two ocean models differ mainly in the representation of the 
vertical coordinate which has important implications for estimat- 
ing mixing, ventilation and the representation of eddies therein. 
Using these two distinct models interchangeably aims to elucidate 
and better assess the sensitivity of the carbon storage in the ocean 
to two numerical representations of the ocean circulation. This 
capability enables us to estimate and address uncertainties in the 
simulated climate due to different ocean physics parameteriza- 
tions. Details in the two ocean models are given below. 

2.2.2. The Russell ocean model 

The Russell ocean model (Russell et al., 1995; Hansen et al., 
2007) is a non-Boussinesq mass-conserving ocean model with 32 
vertical levels and 1 x 1.25° horizontal resolution. All terms are 
solved on an Arakawa C-grid except the Coriolis term that is solved 
on a D-grid. The vertical coordinate is a stretched z-level coordi- 
nate (sometimes described as z*) and has a free surface and natural 
surface boundary fluxes of freshwater and heat. It uses the Gent 
and McWilliams (1990) scheme for isopycnal eddy fluxes and iso- 
pycnal thickness diffusion which is known to lead to weaker and 
more realistic convective mixing, particularly in the Southern 
Ocean. Vertical mixing is parameterized by the K-Profile Parame- 
terization (KPP) vertical mixing scheme (Large et al., 1994) with 
non-local mixing effects included. Vertical advection is discretized 
using a centered differencing scheme. Interactions of ocean with 
other model components as well as tracer advection are done at 
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15 min intervals while the timestep for the model dynamics is 
3.75 min. In this simulation the Bering strait is 3 grid points wide. 

2.2.2. The HYCOM 

The Hybrid Coordinate Ocean Model (HYCOM) is the hybrid ver- 
sion of the Miami Isopycnal Coordinate Ocean Model (MICOM; 
Bleck et al., 1992) and it is described in Sun and Bleck (2006) var- 
iant C. It uses p coordinates near the surface and isopyncals in the 
ocean interior. The version used here has a horizontal resolution of 
1 x 1 cos</>, where </> is the latitude and 26 coordinate layers, a bipo- 
lar grid of matching resolution north of 60°N and an open Bering 
Strait and increased eddy diffusivity. Within 10° north and south 
of the equator, the latitudinal resolution is increased up to 0.03° 
at the equator. 

Each coordinate layer is isopycnal and is assigned a “target” 
density which it follows throughout the ocean except where the 
target isopycnal outcrops. At low latitudes the vertical layers are 
isopycnal in character but become constant-depth layers pole- 
ward. The hybrid coordinate system therefore maintains vertical 
resolution in unstratified water columns, such as in the mixed 


Table 1 

Synopsis of Russell and HYCOM ocean model parameters that are relevant to vertical 
mixing. 



GISSER 

GISSEH 

Atmos-ocean interact, step 

30 min 

30 min 

Tracer adv. timestep 

15 min 

12 h 

Vert. diff. coeff (KPP) 

le-5 m 2 /s 

3e-7 m 2 /s 

Vert. vise, coeff (KPP) 

le-4 m 2 /s 

5e-5 m 2 /s 

Fraction of BL where solar 

1 

Variable, based on 

rad. is absorbed 


Jerlov 

Depth of KPP 

To the bottom of 

To the bottom of 


water column 

mixed layer 

Critical Richardson number 

0.3 

0.45 


layer and also in shallow regions. Judicious choice of target values 
allows the near-surface part of the ocean to always be represented 
by fixed-depth layers. 

A combination of two schemes is used for vertical mixing. In the 
turbulent surface mixed layer, the mixing scheme is KPP (Large 
et al., 1994) with non-local mixing effects included. Beneath the 
mixed layer, a mixing scheme specifically designed for isopycnic 



Fig. 1. Net heat flux (in W m 2 ) climatology and model differences during December-January-February (DJF) and June-July-August (JJA). (a) and (b) observations (OAFLUX 
data), (c) and (d) GISSER and (e) and (f) GISSEH. Positive values mean flux into the ocean. 
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coordinate models (McDougall and Dewar, 1998) is used. In the 
latter, diapycnal diffusivity depends on the smaller of two num- 
bers: a constant diffusivity of 2 x 10 5 m 2 /s and a stratification- 
dependent diffusivity after Gargett (1984). For the tracers, there 
is no isopycnal diffusion. Deep convection at high latitudes is 
parameterized through a convective adjustment scheme. Brine re- 
jected during freezing is uniformly distributed in the uppermost 
200 m if the depth is greater than 1000 m or the total water col- 
umn depth otherwise. The coefficient of the lateral T/S diffusion 
is enhanced near the equator in order to resolve tropical instability 
waves which are known to be highly diffusive. The Gent-McWil- 
liams scheme for parameterizing the effect of baroclinic instability 
is included via direct isopycnal interface smoothing and is limited 
to the isopycnic coordinate domain (which excludes the near-sur- 
face waters). 

The time integration is based on the split-timestep method 
(Bleck and Smith, 1990) in which the barotropic timestep is 50 s 
for the faster modes while the baroclinic time step is half hour. 
The coupling with the atmosphere is synchronous. Tracer transport 
takes place once every 1 2 h. 

Table 1 shows some of the key parameters in the two ocean 
models used here. More details on the physical ocean models can 
be found in Schmidt et al. (2012) as well. 

2.3. The ocean biogeochemistry model, NOBM 

The interactive ocean carbon cycle model consists of a biogeo- 
chemical model (NASA Ocean Biogeochemistry Model, NOBM, 
Gregg and Casey, 2007) and a gas exchange parameterization for 
the computation of the C0 2 flux between the ocean and the 
atmosphere. 

NOBM utilizes ocean temperature and salinity, mixed layer 
depth and the ocean circulation fields and the horizontal advection 
and vertical mixing schemes obtained from the host ocean model 
(Russell or HYCOM) as well as shortwave radiation (direct and dif- 
fuse) and surface wind speed obtained from the atmospheric mod- 
el to produce horizontal and vertical distributions of several 
biogeochemical constituents. 

The abiotic constituents of NOBM (nutrients and carbon) are 
conservative and both biotic and abiotic components are subject 
to different dynamical interactions. There are four phytoplankton 
groups (diatoms, chlorophytes, cyanobacteria, and coccolitho- 
phores), which have distinct maximum growth rates, sinking rates, 
nutrient requirements, and optical properties, in order to represent 
the huge variety of physical and biological environments occurring 
in the global ocean. Four nutrients (nitrate, silicate, ammonium 
and iron), three detrital pools (for nitrates/carbon, silicate and iron) 
and a single herbivore component are also represented. The carbon 
submodel parameterizes the cycling of carbon through the phyto- 
plankton, herbivore and detrital components, affecting the dis- 
solved inorganic and organic carbon in the ocean and interacting 
with the atmosphere. The local profiles of underwater irradiance 
are computed bia a radiative transfer calculation and are needed 
to compute the growth of the phytoplankton groups, their sinking 
profiles due to biologically mediated absorption and scattering in 
the water column. These radiative calculations which depend on 
the spectral nature of the irradiance are based on laboratory exper- 
iments of the optical properties of phytoplankton groups (Gregg 
and Conkright, 2002). 

Particle export in NOBM is represented by gravitational settling 
(sinking). Sinking of particles is variable for the phytoplankton 
groups and detrital components. The corresponding sinking rates 
are derived from mean size estimates for the phytoplankton 
groups, and detritus is adjusted to meet the global distribution of 
nutrients. The values are shown in Table 2. The sinking rate is spec- 
ified for 31 °C and is adjusted for viscosity using Stokes Law (Csan- 


ady, 1986), with temperature parameterization (Gregg and Casey, 
2007; Gregg et al., 2012). 

NOBM has been used previously to understand and evaluate the 
nature of seasonal variability of total chlorophyll, the phytoplank- 
ton functional group and the nutrient distributions in the global 
oceans (Gregg and Casey, 2007) and has been validated extensively 
against in situ and satellite observations (Gregg et al., 2003a). 

2.4. The air-sea gas exchange 

The air-sea gas transfer of C0 2 is parameterized as 

F = k w a([C0 2 ] - [C0 2 ] sot ) (1) 

where F is the flux of C0 2 in/out of the ocean, kw is the gas transfer 
velocity for C0 2 , a is the solubility, a function of pressure, SST and 
SSS, and [C0 2 ] is the surface ocean concentration of C0 2 . [C0 2 ] sat 
is the corresponding saturation concentration of C0 2 in equilibrium 
with a water-vapor-saturated atmosphere at a total atmospheric 
pressure P and a given atmospheric pC0 2 level: 

[C0 2 ] sot = F[C0 2 ]° (2) 

where P° = latm and [C0 2 ]° is the saturation concentration at 
1 atm total pressure. 

The gas transfer velocity is given by, 

/ sc \ 

k "=im) « 2 e> 

where u is the surface wind speed and c is the piston velocity coef- 
ficient taken here equal to 0.337/(3.6 x 10 5 ). The value of c has 
been agreed upon by the Ocean Carbon Model Intercomparison Pro- 
ject, phase II (OCMIP-II) so that the global, annual mean gas transfer 
coefficient for carbon dioxide (k w a) is equal to 0.061 mol/m 2 /yr/ 
patm for preindustrial times. This estimate is based on 
k w apC0 2 = 17 ± 4 mol/m 2 /yr estimated by Broecker et al. (1986) 
based on the global bomb radiocarbon budget divided by the prein- 
dustrial pC0 2 of 280 patm. 5c, the Schmidt number, is computed 
using the temperature of the host ocean model following Wannink- 
hof (1992). k w is computed only over open water. The solubility of 
C0 2 in the water a is also parameterized based on OCMIP using 
prognostic temperature, salinity and sea level pressure. 

In the present set of runs atmospheric concentration of C0 2 is 
held constant at a global average of 285.2 ppmv, although region- 
ally atmospheric C0 2 is allowed to vary due to the distributions of 
the ocean sources and sinks. In addition, air-sea C0 2 is treated as 
radiatively inactive and therefore air-sea ocean C0 2 exchanges 
do not induce any climate change. 

2.5. Radiative coupling of ocean and atmosphere 

In contrast to Gregg and Casey (2007), in the present model, 
NOBM is coupled prognostically to the atmospheric radiation mod- 
ule which is part of modelE2. Direct and diffuse sunlight in 6 wave- 
bands (5 in the near IR and 1 in the visible part of the spectrum) 
computed by the radiation module are decomposed into the 33 
wavelengths (250-3700 nm) needed for the absorption and scat- 
tering by the different phytoplankton species. The decomposition 
is based on the fractions of direct (diffuse) sunlight at the equator 
and each wavelength to the total radiation in each waveband de- 
rived from a monthly climatology provided by the Ocean-Atmo- 
sphere Spectral Irradiance Model (OASIM) and forced by MODIS 
data for aerosols and clouds (Gregg and Casey, 2009). 
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Fig. 2. Surface wind speed (ms 1 ) climatology and model biases with respect to OAFLUX. (a) and (b) DJF and JJA climatologies from OAFLUX, (c) and (d) GISSER bias in DJF and 
JJA, (e) and (f) GISSEH bias in DJF and JJA. 


3. Initialization/evaluation datasets 

3.1. Partial C0 2 pressure in the ocean 

Model estimated pC0 2 in the water, the air-sea flux of C0 2 , sol- 
ubility and piston velocity are compared to the Takahashi atlas 
(Takahashi et al., 2009; Takahashi, 2009) which consists of clima- 
tological mean distributions over the global oceans in non-El 
Nino conditions at 4 x 5 resolution for the reference year 2000. 
The Takahashi atlas is based on pC0 2 measurements obtained from 
1970 to 2007, and utilizes SST and SSS data from the World Ocean 
Circulation Experiment (WOCE) cruises and wind speed data from 
the 1979-2005 NCEP-DOE AMIP-II Reanalysis. Piston velocity is 
parameterized as a function of the square of the wind speed with 
a scaling factor of 0.26 for long-term (monthly) mean gas transfer. 
Wanninkhof DSR-II (2007) suggest a transfer coefficient of about 
0.36 relevant to high-frequency gas exchange. In this work, we 
use a value of 0.26 for the transfer coefficient for consistency with 


previous model evaluation studies based on OCMIP-II parameter- 
izations although perhaps Wanninkhof s latest estimate would be 
more appropriate in a coupled model context where coupling is 
synchronous and shorter rather than longer time scales affect the 
flux. The uncertainties in the Takahashi atlas are estimated at 
30% and are higher than the long term mean. 

Data from the study areas/stations such as the Bermuda Atlantic 
Timeseries Study (BATS, ftp://ftp.bios.edu/BATS), the Hawaii Ocean 
Timeseries Study (HOTS, Dore et al., 2009) and CARIOCA (http:// 
www.imars.usf.edu/CAR/index.html) are also used for the valida- 
tion of the seasonal cycle of pC0 2 in the models. 

We note that, the comparison with the Takahashi atlas is used 
here only qualitatively to depict the spatial patterns of pC0 2 in 
the ocean. For the purposes of this paper, the Takahashi atlas is 
not an exact validation tool, since the model runs described here 
refer to a preindustrial atmosphere, in contrast to the Takahashi at- 
las which describes the state of ocean surface pC0 2 for atmospheric 
concentrations in year 2000. 
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3.2. DIC 

Initial distributions of the dissolved inorganic carbon are taken 
from the Global Ocean Data Analysis Project (GLODAP). These dis- 
tributions are the result of a synthesis of approximately 8 yr of data 
(1990-1998) from the global C0 2 survey that included 50 individ- 
ual cruises. The data were unified into an internally consistent 
dataset (Key et al., 2004). 

In particular, for the preindustrial spinup and the equilibrium 
runs described here, estimates of the preindustrial (i.e. present 
day - anthropogenic) distributions of DIC were used. Deviations 
of model results at equilibrium from the GLODAP data are useful 
in explaining the different way processes are resolved in the two 
models. 

3.3. Nutrients, SST, SSS, mixed layer depths 

Nutrients (nitrates and silicates) initial conditions, and SST and 
SSS evaluation datasets are obtained from the National Oceano- 


graphic Data Center (NODC) World Ocean Atlas 2005 (Locarnini 
et al., 2006; Antonov et al., 2006; Garcia et al., 2006). These data 
are provided in a monthly mean climatology at 1° resolution and 
33 nominal depths and are interpolated to each ocean model’s hor- 
izontal and vertical grid. 

Model mixed layer depths are evaluated against a 2° global cli- 
matology which is constructed using individual T, S profiles and a 
combination of the critical difference values for temperature and 
density from the surface value, 3T = 0.2 °C or Sp = 0.03 kg m“ 3 
(deBoyer Montegut et al., 2004). 

3.4. Dust-iron 

Following Gregg et al. (2003b), the dust-iron deposition flux is 
obtained from the Georgia Institute of Technology/Goddard Global 
Ozone Chemistry Aerosol Radiation and Transport climatology for 
the period 1982-2000, except for the 1997-1999 period when data 
was unavailable (GOCART; Ginoux et al., 2001). The iron content is 
assumed to vary within each dust size class, with clay containing 



Fig. 3. Sea surface temperature climatologies (°C) from NODC and model biases, (a) DJF SST from NODC, (b) JJA SST from NODC, (c) and (d) GISSER biases from NODC during 
DJF and JJA respectively and (e) and (f) GISSEH biases from NODC during DJF and JJA respectively. 
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3.5% iron and silt containing 1.2% iron (Fung et al., 2000). Further- 
more, a constant iron solubility of 2% is also assumed. 

3.5. Surface radiative and thermal fluxes and surface wind speed 

Surface radiative fluxes (shortwave and longwave) are obtained 
from the International Satellite Cloud Climatology Project (ISCCP) 
FD flux dataset. The resolution is 3 h at 2.5° covering July 1983 
through December 2007 (see Zhang et al., 2004 for details and ref- 
erences). Uncertainties of the monthly fluxes are estimated to be in 
the range of 10-15 W m -2 at the surface. 

Latent and sensible heat fluxes come from the Woods Hole 
Oceanographic Institute (WHOI) Objectively Analyzed air-sea Heat 
Fluxes (OAFlux) surface heat flux dataset and are used here for 
model evaluation. The daily data cover the period 1985 to 2009, 
while the monthly fluxes are available from 1958 onward. OAFlux 
is a blended dataset created by combining satellite-derived obser- 
vations of meteorological variables with those from numerical 
weather prediction (NWP) model output. Comparison to in situ 
flux measurements shows that the mean average of the differences 


between OAFlux and buoys over 107 locations (a measure of bias) 
is l.OWm" 2 , and that the absolute mean average of the OAFlux- 
buoy differences (a measure of variance) is 7.4 W m -2 . For further 
information, see Yu et al. (2008). 

4. Results 

In the results presented below, all runs start from equilibrium 
conditions obtained after about 3000 yr of integration of the cou- 
pled atmosphere-ocean-ice models. The upper ocean is in statisti- 
cal steady state with respect to all thermodynamic variables 
whereas the deeper ocean is still evolving, albeit very slowly and 
on timescales longer than the ones of interest here. Both models 
are then coupled to NOBM, using initial conditions for nutrients 
from NODC and Fung et al. (2000) and DIC from the GLODAP pre- 
industrial estimate (Key et al., 2004). Regional variations due to the 
ocean-atmosphere flux are computed at each grid point but the 
global average atmospheric C0 2 concentration is held constant at 
285.2 ppmv. The models are further integrated until the C0 2 flux 
at the surface of the ocean is in statistical steady state and approx- 


(a) Potential Temperature 15N-90N 


(b) Potential Temperature 15S-15N 




(c) Potential Temperature 15S-90S 



Fig. 4. Potential temperature profiles (°C) averaged over latitude zones, from the two models and the NODC data, (a) 15-90° N, (b) 15° S-15° N and (c) 15-90° S. 
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imately 0.9PgCyr _1 . This corresponds to an additional 300 yr of 
integration after which each coupled ocean-atmosphere-biogeo- 
chemistry model is in equilibrium. 

In the following discussion, observations are used only as a 
benchmark to assess the magnitude and realism of the intermodel 
differences, since observations refer to present day climate condi- 
tions whereas models simulate here an unperturbed, stable, prein- 
dustrial climate. 

4.1. Physical ocean model differences 

Surface heating (i.e. the sum of radiative and thermal surface 
fluxes into the ocean) is shown in Fig. 1. Both model estimates 
agree well with the observed seasonal cycle and the inter-hemi- 
spheric differences but exaggerate the seasonal cycle amplitude 
(Fig. 1). In the subtropics and midlatitudes of each hemisphere, 
the models overestimate heat losses during the cold season (e.g. 
December-January-February (DJF) for the Northern Hemisphere 
and June-July-August (JJA) for the Southern Hemisphere) and 


underestimate the heat uptake during the warm seasons. In the 
tropics, both models show weaker than observed heat gain during 
both seasons (Fig. 1), with the eastern Equatorial Pacific and the 
tropical Atlantic receiving about 50-100 W m -2 less than what OA- 
FLUX shows. The Pacific cold tongue region is weaker and narrower 
in both models compared to observations. 

The model differences in surface heating patterns are related to 
clouds, sea-ice thickness and extent as well as the SST distribu- 
tions. Both models overestimate cloud amount over the tropical re- 
gions by about 30-50% (Schmidt et al., 2012) resulting in reduced 
heat reaching the surface there. The models also overestimate the 
wintertime (DJF) cloud amount over the Arctic Ocean (Schmidt 
et al., 2012) thereby limiting the amount of surface cooling. In 
the Southern Ocean, though, during the austral winter (JJA) both 
models underpredict cloud cover, which leads to extreme (about 
200 W m -2 ) heat losses near the coast of Antarctica. 

Surface winds are weaker in the GISS models by about 1- 
3 m s -1 (Fig. 2(a)-(d)) compared to the OAFLUX data, with the lar- 
ger differences exhibited by the GISSER model. The tropical regions 



Fig. 5. Surface salinity climatology from NODC and model differences, (a) and (b) SSS from NODC for DJF and JJA respectively, (c) and (d) SSS biases in GISSER for DJF and JJA 
respectively and (e) and (f) SSS biases in GISSEH for DJF and JJA respectively. 
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as well as the wintertime storm tracks (in the northwest Atlantic 
and Pacific Oceans and in the Southern Oceans) are the regions 
where the largest underestimates exist. In contrast, both models 
overpredict the coastal winds along the western continental bor- 
ders during both DJF and JJA, due to the inability of the low resolu- 
tion grid to resolve steep topography there. 

As a result, surface temperature is about 2-4 °C warmer than 
observations in both models in the Southern Ocean, with the larger 
overestimates displayed by GISSEH during the warm season (DJF, 
Fig. 3). These differences are explained by the lack of cloudiness 
that allows more heating to reach the surface, increasing the ther- 
mal stratification near the surface, further inhibiting vertical ex- 
changes which would otherwise cool the surface. In the Northern 
Hemisphere, the colder wintertime SSTs are linked to increased 
heat losses due to the decreased cloudiness. In the tropics, the 
much warmer surface ocean in the models, which is due to ocean 
model deficiencies, is compounded by the increased clouds which 
prevent radiative cooling of the ocean surface. 

Zonally averaged potential temperature vertical profiles (Fig. 4) 
reveal how model differences are transferred to the deeper ocean 


levels. In the 3 regions depicted (Northern Hemisphere extratrop- 
ics, tropics and Southern Hemisphere extratropics), model esti- 
mates bracket the observations down to about 200 m, with 
GISSER being colder than the observations and GISSEH. Between 
200 m and 3500 m depth, GISSER temperatures are much warmer 
than both observations and HYCOM (Fig. 4). Such mid-depth 
warming in GISSER is a result of the excessive warming of the sur- 
face ocean in the deep convection regions, which forms waters that 
ventilate the mid-ocean depths with warm biases (seen in Fig. 3). 

Surface salinity differs significantly between the two models 
and from observations particularly in the Arctic Ocean (Fig. 5). 
Large overestimates exist there because of increased brine rejec- 
tion during excessive ice formation in DJF as well as due to defi- 
cient freshwater river outflow (Schmidt et al., 2012). In the 
summer months (JJA) in the Arctic Ocean both models underpre- 
dict the sea ice coverage, which implies excessive melting of the 
wintertime sea ice leading to lower surface salinity. In the South- 
ern Ocean, GISSER overestimates sea surface salinity compared to 
GISSEH, because of the larger evaporation in that model (Fig. 6). 
Surface salinity is also overestimated in the models at the coastal 


GISSER: Evaporation, JJA 
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Fig. 6. Evaporation (in W m 2 ) climatologies from observations and the two model runs for DJF (left columns) and JJA (right columns). 
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upwelling regions again due to excessive evaporation in these re- 
gions, more clouds and warmer SSTs. Weaker than observed coast- 
al winds impose limited upwelling, drawing to the surface warmer 
water from shallower depths, thereby enhancing evaporation and 
leading to increased SSS distributions. 

Elsewhere in the tropics and in the subtropics, slightly lower 
than observed model surface salinities are attributed to the in- 
creased model precipitation compared to GPCP climatology 
(Schmidt et al, 2012) which exceeds observations in the tropics 
by about 2-10 mm d -1 and the subtropics by 1-2 mm d _1 . 

Mixed layer depths (MLDs) are shallower than observed in the 
Southern Ocean in both models, especially in GISSEH (Fig. 7). Deep 
convection in the North Atlantic reaches greater depths in GISSEH 
but is confined in a smaller geographical region in the models than 
in the observations, leading to a decreased mid-depth warming 
bias in GISSEH (Fig. 4), and therefore more realistic deep ocean 
temperature profiles than GISSER. During JJA, wintertime mixed 
layer depths against the Antarctic coast are deeper than observed 
in GISSER due to reduced cloudiness and too little sea-ice, exces- 
sive cooling and dense water formation. By contrast, in the Arctic 


Ocean, winter mixed layer depths in both models are shallower 
than observed due to the insulating effect of too much model ice. 

4.2. The air-sea surface C0 2 flux: solubility pump differences 

During spinup, both ocean models emit C0 2 to the atmosphere 
for an adjustment period of about 250 yr, implying that the oceanic 
initial condition has too much surface DIC and is out-of-balance 
with the model atmosphere. The rearrangement of DIC within 
the mixed layer, together with the slow evolution of deep ocean 
properties are responsible for this adjustment period, since sea sur- 
face temperature, surface winds or mixed layer depth are in statis- 
tical steady state while the air-sea C0 2 flux is not interactive with 
the model’s radiation and does not affect the model climate. 

The spatial distribution of the C0 2 flux in the two models at 
equilibrium for a climatological December-January-February 
(DJF) and June-July-August (JJA) is shown in Fig. 8, together with 
the distribution of the flux from the Takahashi atlas. Both models 
predict a strong sink of C0 2 in the cold season hemisphere. During 
DJF, the North Atlantic deep convection region (Labrador Sea, Ice- 


(a) DeBoyer: Mixed Layer Depth, DJF 


(b) DeBoyer: Mixed Layer Depth, JJA 
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Fig. 7. Mixed layer depths climatologies (in m) from (a) and (b) the deBoyer atlas, (c) and (d) GISSER and (e) and (f) GISSEH for DJF and JJA. 
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land and Norwegian basins) dominates the global C0 2 drawdown 
(Fig. 8(c) and (e)) which is more pronounced in GISSEH than in GIS- 
SER. The maximum drawdown in the N. Atlantic sink region is 
about 8 mol C0 2 m -2 yr -1 in GISSER and about 6 mol C0 2 m -2 yr -1 
in GISSEH although the latter over a more extended area. As shown 
in the previous section, DJF SSTs are warmer by about 5 °C in GIS- 
SEH than GISSER (Fig. 8(c) and (e)) therefore SST does not drive the 
more vigorous uptake in GISSEH. Storm tracks, on the other hand, 
are stronger in GISSEH than in GISSER (Fig. 2(c) and (e)) by about 
1-2 m s -1 on the monthly mean which partly explains the stronger 
C0 2 sink in this model. At the same time, surface ocean pC0 2 in 
GISSEH is lower (Fig. 9) which further enhances the C0 2 uptake 
by the GISSEH ocean moreso than the GISSER ocean. It is therefore 
the dependency of the flux (Eqs. (l)-(3)) on the surface pC0 2 and 
the surface wind speeds that mainly determines the uptake in 
the North Atlantic rather than the SST field in these two models. 

During JJA, the North Atlantic remains a sink region of C0 2 , 
although a little weakened in magnitude because the storm tracks 
are weaker (Fig. 8(d) and (f)). Additionally, the ocean between 20°S 


and 50°S becomes a major carbon sink (Fig. 8(d) and (f)) with a 
max of about 2 mol C0 2 m -2 yr -1 in GISSER and about 6 mol C0 2 - 
m -2 yr -i j n GISSEH. This elevated uptake in the GISSEH Southern 
Ocean compared to GISSER relates to the higher surface pC0 2 in 
GISSEH (by about 50 ppmv) than in GISSER (Fig. 9(c) and (d)) de- 
spite the warmer SSTs and the weaker surface winds in GISSEH. 

At the equator, the Eastern Pacific source of C0 2 is more pro- 
nounced in the GISSER model than the GISSEH during both seasons 
(Fig. 8) due to greater pC0 2 and stronger trade winds (Fig. 2). Both 
these effects overpower the significantly warmer SST in the GISSER 
model, particularly along the equator. Elsewhere in the tropics, 
such as in the Warm Pool region, GISSEH pC0 2 is much higher than 
the atmospheric value (285.2 ppmv; Fig. 9) and therefore the out- 
gassing is more extensive. Here, the very high values of pC0 2 are 
related to the very low salinities (Fig. 5). 

Annually, the Southern Hemisphere (south of 20°S) emits about 
0.11 PgCyr- 1 in GISSER and O.ayPgCyi- 1 in the GISSEH. Compared 
to the Takahashi atlas (Takahashi et al., 2009) in which the same 
region uptakes about 0.08 PgCyr -1 , the larger and opposite sign 



Fig. 8. Flux of C0 2 into the ocean (mol-C0 2 m 2 yr 1 ) in the (a) and (b) Takahashi 2009 atlas, (c) and (d) GISSER and (e) and (f) GISSEH models for DJF (left panels) and JJA 
(right panels). Positive flux denotes ocean uptake. 
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flux in the GISS models is attributed to the larger disequilibrium of 
pC0 2 in the preindustrial climate that is modeled here than in the 
present climate (and hence the Takahashi atlas). In the equatorial 
regions GISS models predict a source of about 1.31 PgCyr -1 in 
the GISSER and 1.01 Pg Cyr -1 in the GISSEH. The equatorial source 
in the Takahashi atlas is 0.69 Pg Cyr -1 again due to the lower atmo- 
spheric concentrations of C0 2 in the model climate. The Northern 
Hemisphere oceans uptake about 0.45 Pg Cyr -1 in the GISSER and 


0.46 Pg Cyr -1 in the GISSEH annually, whereas in Takahashi the up- 
take is estimated at 1.3 Pg Cyr -1 . 

The geographical patterns of the natural C0 2 flux into the ocean 
compare well, both qualitatively and quantitatively, with inverse 
estimates from about 10 natural carbon cycle models as described 
in Mikaloff-Fletcher et al. (2007). Fig. 10 shows the GISSER and GIS- 
SEH flux integrals over several geographic regions as well as the 
multi-model mean and standard deviation from the inverse esti- 
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Fig. 10. C0 2 flux (in Pg Cyr 1 ) into the ocean integrated over different regions, as in Mikaloff-Fletcher et al. (2007). 
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pC0 2 anomaly at BATS pC0 2 anomaly at HOTS 



Fig. 11. pC0 2 values from observations, the Takahashi atlas and GISSER and GISSEH models at the locations of the Bermuda Atlantic Timeseries Study area (BATS) and the 
Hawaii Ocean Timeseries Study (HOTS). 





Fig. 12. Longitudinal variation of the depth integrated DIC tendencies (mmol Cm -2 h -1 ) as described by different terms in Eqn. 5 in GISSER only, along lines of constant 
latitude (top) 50° N, (middle) the equator, and (bottom) 70° S. Solid blue lines represent phytoplankton respiration, dashed blue lines represent remineralization, cyan lines 
bacterial degradation of DOC, green lines zooplankton respiration, solid red lines show gas exchange, and dashed red lines show phytoplankton growth. Dashed black lines are 
the sum of these biology terms and solid black lines show advection (mean + eddy). Note vertical axes range is different in the three panels. 


mates. Both GISS models predict the natural C0 2 flux within the 
uncertainty range of the multi-model mean in most regions of 
the globe. Most notably, both GISS models fall outside the uncer- 
tainty ranges in the polar Southern Ocean where as Mikaloff- 
Fletcher et al. (2007) discuss even inverse model estimates are 
highly uncertain, and in the Arctic Ocean and the North Atlantic 
tropics. 

In situ data in stations HOTS and BATS provide a useful tool 
assessing the quality of the modeled seasonal cycle of C0 2 . Signif- 
icant difficulty arises, however, when one tries to compare mod- 
eled fields which are smooth functions over a gridsize area (here 
about 100 km) with point measurements. Nevertheless, Fig. 11 
shows that both GISS models capture the phase and magnitude 


of the seasonal cycle exhibited in the Takahashi atlas, but underes- 
timate its magnitude significantly at these locations. 

Surface ocean pC0 2 biases reflect the DIC differences over most 
regions which will be discussed next. 

4.3. The air-sea C0 2 flux: biological pump differences 

DIC changes in the ocean are due to the air-sea flux of C0 2 , met- 
abolic processes, vertical mixing and lateral eddy diffusion (Gregg 
and Casey, 2007): 

f t [D/C] = V • (KV[D/q) - V • V[D/q - <DX,//,P, + mZiUfi (4) 

+ <D©H + ^[DOq+Pa N D c /(C:N)+F C o 2 (5) 
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Table 2 

Vector sinking rate of phytoplankton and carbon detritus 
at 31 °C, in m d _1 (Gregg et al., 2012). 


Sinking rate (m d 1 ) 

Diatoms 

0.75 

Chlorophytes 

0.25 

Cyanobacteria 

0.0085 

Coccolithophores 

0.3-1. 4 

Carbon detritus 

40 


Table 3 

Regional correlations between modeled and observed annual climatologies of total 
chlorophyll (tchl), diatoms (diat), chlorophytes (chlo), coccolithophores (cocc), nitrate 
(nitr), silicate (sili), iron and DIC in the two models. 



tchl 

diat 

chlo 

cyan 

cocc 

nitr 

sili 

iron 

die 

GISSER 

0.64 

0.78 

0.86 

0.77 

0.39 

0.98 

0.36 

0.25 

0.84 

GISSEH 

0.52 

0.75 

0.86 

0.59 

0.74 

0.97 

0.61 

0.18 

0.50 


where the first term on the right-hand-side is the diffusion of DIC 
(turbulent mixing) in which I< is parameterized based on Large 
et al., 1994. The second term is horizontal advection which includes 
eddy transport, followed in turn by phytoplankton (P,-) growth and 
respiration, zooplankton (H) respiration, bacterial degradation of 


dissolved organic carbon (DOC), remineralization of detrital carbon 
and finally, F C02 , the flux due to gas exchange between the ocean 
and the atmosphere. 

Fig. 12 (top panel) shows that gas exchange (solid red line) 
along 50°N is the leading term in the DIC changes in the North 
Atlantic deep convection region (longitudes: -70E to -40E) where 
intense C0 2 uptake occurs. Everywhere else at that latitude, chlo- 
rophyll growth (dashed red line) dominates all other terms in the 
DIC equation. Along the equator, phytoplankton growth (dashed 
red line) becomes the main driver for DIC decreases, and is further 
reinforced by outgassing. Similarly in the Southern Ocean, near the 
Antarctic coast (70°S), phytoplankton growth is the main sink of 
DIC. Horizontal advection (sum of the mean advection and the 
eddy stirring term) is of the same order of magnitude, opposite 
sign as the total biology terms nearly balancing the local DIC 
changes as expected from Eq. (5). 

The biologically mediated DIC uptake is larger in GISSER be- 
cause total chlorophyll is higher in this model over most regions 
except in the Southern Ocean south of about 50°S (Fig. 13). The 
geographical pattern is similar during both DJF and JJA although 
the differences are more pronounced during the warm season 
(DJF in the Southern Ocean; JJA in the Arctic). 

During DJF, in North Atlantic higher chlorophyll distributions in 
GISSER (Fig. 13(a)) are sustained by increased nitrate load 
(Fig. 14(a)) which is due to deeper mixing there (Fig. 7), despite 
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the significantly colder SSTs (Fig. 3). On the contrary, the JJA 
blooms in the Arctic (Fig. 13(b)) are related to the increased warm- 
ing of the water column (Fig. 3) as more ice melts in GISSER than in 
GISSEH. 

In the tropical regions, over a narrow latitude band around the 
equator, GISSER has lower chlorophyl concentrations during both 
seasons, due to the colder SSTs as well as lower concentrations of 
nitrates (Fig. 14(a) and (b)). 

In the Southern Ocean, south of 50°S particularly during DJF 
(the austral summer) chlorophyll is lower in GISSER because of 
lower nitrate load (Fig. 14(a)) in this model. Shallower mixed lay- 
ers (Fig. 7) result mixing that does not extend deep enough to bring 
more nutrients to the surface. 

In the Antarctic convergence zone (around 45°-50°S) increased 
biological production in GISSER is associated to the abundant ni- 
trates there, due again to deeper MLD even though SST is colder. 
This is the area, the Subantarctic Front, where the Subantarctic 
Mode waters form, that is characterized by high nitrate and low 
silicate load and increased productivity. Surface waters in this re- 
gion originate from the deep North Atlantic and are rich in remin- 
eralized nutrients. 


Correlations between observed and modeled climatologies in 
each region (Table 3) reveal that GISSER has better skill in simulat- 
ing ocean productivity and main phytoplankton groups than GIS- 
SEH. GISSER, and less so GISSEH, captures total chlorophyll and 
nitrate spatial patterns with high fidelity however, coccolitho- 
phores and silicate/iron distributions are not very realistic. 

Deep carbon storage and export can be assessed through the 
model DIC biases. The annual-mean, zonal distributions of DIC in 
the two models are shown in Fig. 15. Surface DIC is lower in the 
equatorial regions in GISSER because of the more vigorous outgas- 
sing (Fig. 8). In the Northern Hemisphere subpolar regions, again, 
surface DIC is lower in the GISSER model because of biological up- 
take since there are extensive chlorophyll blooms there throughout 
the year (Fig. 13). In the Southern Ocean, however, surface DIC is 
higher in GISSER due to more intense mixing with deeper waters 
that are richer in DIC (Fig. 15) whereas the GISSEH sharp carbocline 
prevents such mixing. As a result, deep DIC distributions are higher 
in the GISSEH model (Fig. 16). 

Comparison of the deep DIC distributions in the two models 
with the initial conditions (GLODAP) (Fig. 16) shows that at 
1500 m GISSEH holds more DIC than GISSER and both models have 
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Fig. 15. Zonally averaged annual mean DIC in mmol Cm 3 . On the left are distributions from GISSER and on the right from GISSEH. 



Fig. 16. DIC differences between the models and the GLODAP estimates at 1500 m. (a) GISSER-GLODAP and (b) GISSEH-GLODAP. Units are in mmol Cm 3 . 
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Fig. 18. Alkalinity climatology and model differences, (a) Alkalinity from GLODAP, (b) GISSER-GLODAP and (c) GISSEH-GLODAP. 


more DIC than GLODAP particularly in low latitudes. Those pat- 
terns are explained by the more vigorous uptake in GISSEH and 
the higher productivity in GISSER in the North Atlantic convection 
region. Both models show depletion of DIC in the deep Southern 
Ocean regions, due to the intense outgassing and the deep mixing 
in these regions in both models, especially in HYCOM. 

In addition to ventilation from the surface, another source for 
DIC at depth is the remineralization of the carbon detritus. GISSER 
has more detritus over most extratropics (Fig. 17) because of colder 
temperatures (Fig. 4) and less efficient remineralization there (Eq. 
(5)). Conversely, at low latitudes GISSER is warmer which implies 


increased remineralization rates and leads to lower detritus con- 
centrations. As a result, the globally integrated carbon export at 
75 m is about 6.7 Pg Cyr -1 in GISSER and about 2 Pg Cyr -1 in GIS- 
SEH. In addition to weaker remineralization in GISSEH, other rea- 
sons for the decreased carbon export in GISSEH are the lower 
chlorophyl concentrations, the decreased vertical mixing and the 
subsequent lower surface DIC concentrations, as well as the more 
intense outgassing. 

Since primary productivity is about 17 PgCyr -1 in GISSER and 
2 Pg Cyr -1 in GISSEH, the carbon export efficiency, defined as the 
ratio of the globally integrated export and the globally integrated 
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primary productivity are very similar; 0.4 for GISSER and 0.3 for 
HYCOM. 

4.4. The air-sea surface C0 2 flux: carbonate pump differences 

The air-sea flux differences between the two models are 
strongly correlated with the pC0 2 differences, as discussed in Sec- 
tion 4.1. In addition to DIC, SST, SSS and surface nutrient distribu- 
tions discussed previously, alkalinity also affects the DIC budget via 
the air-sea gas exchange. Theoretically, alkalinity changes with 
precipitation and evaporation as well as the production and disso- 
lution of calcium carbonate, CaC0 3 , and the assimilation and re- 
lease of dissolved inorganic nitrogen. However, in the present 
model formulation, alkalinity is a constant and uniform function 
of salinity, and does not affect pC0 2 except only at the surface 
ocean where it is involved in gas exchange. Therefore, salinity 
biases discussed before (Fig. 5) reflect large underestimates of alka- 
linity in the Southern Ocean mainly (Fig. 18), which are more pro- 
nounced in the GISSEH model. In the North Atlantic, alkalinity 
biases are such that they explain the lower C0 2 uptake in the GIS- 
SER model during DJF whereas in the Southern Ocean alkalinity 
biases lead to more outgassing in the GISSEH. Compared to GLO- 
DAP estimates, model alkalinity is overpredicted in the upwelling 
coastal regions but is underestimated everywhere else, particularly 
in the Southern Ocean. This implies that in this region the exces- 
sive model outgassing compared to Takahashi (Fig. 8), cannot be 
explained by the alkalinity biases there. 

5. Conclusions 

This paper documents the similarities, differences and possible 
driving mechanisms of the natural carbon cycle between the two 
NASA-GISS climate models. These two climate models share the 
same atmospheric, land, ice and ocean carbon components but dif- 
fer in the ocean physics component. The ocean models employ dif- 
ferent vertical coordinate system which is known to produce 
different distributions of seawater properties at depths where 
numerical and spurious diapycnal mixing can become important 
(Bleck et al., 1992; Griffies et al., 2000). Moreover, the ocean mod- 
els differ in the details of the vertical mixing scheme (KPP); in GIS- 
SER KPP extends to the bottom of the water column but non-local 
terms are not included whereas in GISSEH it only reaches the bot- 
tom of the mixed layer. Eddy mixing along isopycnals (Gent and 
McWilliams, 1990) is only explicitly modeled in GISSER, whereas 
in GISSEH is approximated by isopycnal thickness diffusion for 
the regions where these isopycnals do not outcrop. Moreover, 
there is no tracer horizontal advection in GISSEH. Although the 
models are different in other ways, the significantly decreased ver- 
tical diffusion/mixing below the mixed layer in the GISSEH is per- 
haps the main difference between them. 

Model results, even though they refer to an unperturbed, stable 
preindustrial physical climate, are assessed against observations 
from present day climate, in order to estimate how large inter- 
model differences are compared to the anthropogenic signal, there- 
by setting the stage for present-day GISS model evaluations. 

We find that model differences in the natural C0 2 flux are ex- 
plained by differences in surface temperature, salinities and winds, 
which are in turn controlled by the different distributions of clouds 
and ice cover in the two models as well as processes in the interior 
of the ocean. We conclude that the air-sea flux is most sensitive to 
differences in the strength of the surface winds in the two models 
and the differences in the partial C0 2 distributions, in general 
agreement with Doney et al. (2004, 2006). 

Storm-track winds in the North Atlantic are stronger in the GIS- 
SEH model, due to stronger evaporative feedbacks driven by the 


warmer SSTs. Lower vertical, interior ocean mixing in GISSEH, 
though, which is an inherent feature of layer models, overwhelms 
the wind effect, slowing down the thermohaline circulation in the 
model, resulting in global northern hemisphere meridional over- 
turning circulation (MOC) of about 18 Sv as opposed to 25 Sv in 
GISSER. Deep convection in GISSEH reaches shallower depths than 
GISSER and therefore brings to the surface lower amounts of DIC 
and nutrients. Total chlorophyll and primary productivity are 
therefore larger in GISSER in the North Atlantic due to more nutri- 
ents. Therefore the more intense uptake in HYCOM is solubility dri- 
ven (winds, pC0 2 ) rather than biologically driven (primary 
production). 

In the tropics, stronger outgassing in the GISSER model, partic- 
ularly in the eastern equatorial Pacific, is attributed to higher sur- 
face temperatures but also stronger trade winds in this model. Still, 
equatorial upwelling and outgassing is reduced in GISSER com- 
pared to observations (Takahashi atlas) due to weaker than ob- 
served trade winds. This is a complex model deficiency 
compounded by the lack of horizontal resolution in the Russell 
ocean model which therefore does not resolve the equatorial 
waveguide and the associated baroclinic instabilities. However, 
other differences between the ocean models, such as horizontal 
advection and eddy stirring are important in this region and are 
modeled differently in the two models. 

The Southern Ocean is a region where the large scale winds 
through Ekman layer dynamics subduct most of the carbon dioxide 
that enters the ocean (Ito et al., 2010) while eddies, which are 
responsible for the lateral transport through eddy stirring, play a 
lesser role. We find that the models’ Southern Ocean fluxes are of 
considerably different magnitude and/or sign from the Takahashi 
atlas, and that these differences cannot be explained by alkalinity, 
SSTs or surface winds. Both models cannot resolve adequately lat- 
eral eddy mixing resulting in mixed layer depths with are too shal- 
low in this region, especially in GISSEH, leading to unrealistic 
subduction processes. Moreover, the particularly shallow mixed 
layer depths in HYCOM are a result of the higher SSTs and the po- 
sitive SST-MLD feedback in which shallower mixed layers lead to 
more heating of the surface layer, thereby increasing the static sta- 
bility of the water column and further reducing vertical mixing. 
The higher SSTs affect gas exchange not through solubility but 
rather through the shallow and sharper thermocline and the re- 
duced exchanges of DIC and nutrients with the deeper ocean. 

In the Southern Ocean two distinct regions emerge; north of 
50°S, i.e. roughly north of the divergence zone, where the GISSEH 
ocean surface is warmer and slightly more saline, winds are weak- 
er (especially during DJF), pC0 2 is greater, total chlorophyll and 
nutrients are lower particularly along the Subantarctic Conver- 
gence Zone and south of 50°S, where the GISSER ocean model outg- 
asses less, is warmer and more saline at the surface, winds are 
weaker total chlorophyll and nutrients are lower. In the GISSER 
model the Antarctic Bottom Water (ABW) forms near the Antarctic 
coast and unrealistic mixing processes at depth in GISSER result in 
the fast depletion of ABW. The area between SAZ and the diver- 
gence zone is very sensitive to parameterizations of vertical mixing 
and as well as the isopycnal diffusion mixing (Gent-McWilliams) 
and since these are modeled differently in our two ocean models, 
this region emerges as the most sensitive to ocean model parame- 
terizations. This is in agreement with the idea that excessive 
Southern Ocean mixing (vertical diffusion or eddy stirring) and/ 
or weaker winds (as they are in both models here, especially GIS- 
SER), the dominant transformation of abyssal waters into light 
water takes place within the low-latitude pycnocline instead of 
in the Southern Ocean (Gnanadesikan et al., 2002). Such a regional 
difference greatly affects the Southern Ocean sink of C0 2 in both 
models, particularly in GISSER. Nevertheless, the Southern Ocean 
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region is one of highly uncertain estimates of the C0 2 gas exchange 
in several other models, as shown in Mikaloff-Fletcher et al. (2007). 

GISSER has better skill in simulating the biological pump than 
GISSEH. Very little mixing between the upper ocean and the inte- 
rior brings less nutrients to the surface reducing the biological pro- 
duction in GISSEH. GISSER, and less so GISSEH, captures total 
chlorophyll and nitrate spatial patterns with high fidelity however, 
coccolithophores and silicate/iron distributions are not very realis- 
tic. The lack of processes such as denitrification and the highly 
approximate alkalinity balance in both models are probably 
responsible for these biases. 

The two models tend to have the largest DIC values at mid- 
depths, along the NADW pathway. Therefore, GISSEH which is 
characterized by reduced interior mixing, is colder and holds more 
DIC there at these depths. In the Southern Ocean, both models have 
lower DIC than GLODAP due to deeper mixing and large outgassing 
there. 

Both models have lower than observed primary productivity 
and carbon export, despite the fact that they simulate fairly well 
total chlorophyl concentrations, which is a shared deficiency 
among other models (Carr et al., 2006). We are presently investi- 
gating the sensitivity of the carbon cycle model (NOBM) to remin- 
eralization and sinking parameterizations. GISSEH exhibits 
particularly reduced primary productivity and export compared 
to GISSER both of which result from a combination of factors such 
as the lower chlorophyl concentrations, the decreased vertical mix- 
ing and the subsequently reduced surface DIC and nutrients con- 
centrations. Carbon export efficiency though, which is a measure 
of the strength of the biological pump, is similar in the two models, 
since the same processes which are responsible for keeping pri- 
mary productivity low are also responsible for keeping carbon ex- 
port low. 

To conclude, physical ocean model differences which feed into 
atmospheric fields, such as the storm track strength, evaporation 
and precipitation and net heating over the ocean, result in solubil- 
ity and hence uptake differences in the two models. Vertical mix- 
ing differences are mostly responsible for the divergent model 
responses with regards to the biological pump. Therefore the qual- 
ity of the ocean component is critical for the realistic representa- 
tion of the natural C0 2 ocean sink. Future studies will assess how 
the differences between the ocean model components might im- 
pact the sensitivity or the transient behavior of the oceanic carbon 
cycle in the GISS climate models. 
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